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MATHEMATICAL  ANALYSIS  OF  ARTIFICIAL 
ENZYMATICALLY  ACTIVE  MEMBRANES 

Preface 

This  report  contains  material  from  lectures  given  in  May 
1979  at  Brown  University.  The  systems  described  here  are  im¬ 
mobilized  enzyme  membranes,  produced  and  experimentally  studied 
by  D.  Thomas  and  his  group  (Laboratoire  de  Technologie  Enzyma- 
tique,  E.R.A.338,  Universite  de  Technologie  de  CompiSgne,  Compifcgne, 
France) .  The  concentrations  in  these  artificial  membranes  are 
governed  by  the  interaction  of  diffusion  and  reaction,  and  they 
may  have  behaviors  characteristic  of  dissipative  structures: 
hysteresis  (related  to  the  existence  of  multiple  stable  steady 
states),  oscillations,  wave  front  propagation  and  pattern  forma¬ 
tion.  Since  analogous  phenomena  exist  in  (much  more  complex!) 
living  systems  (short-term  memory, biological  clocks,  transmission 
of  signals,  morphogenesis)  it  is  interesting  to  study,  in  a  well- 
defined  context,  the  interaction  of  diffusion  and  enzyme  reaction. 
Further  details  and  more  references  on  diffusion-reaction  phenomena 
can  be  found  in  the  books  of  Aris  [1],  Nicolis  and  Prigogine  [2], 
Murray  [3]  and  Fife  [4]. 

The  author  wishes  to  express  his  thanks  to  the  members  of 
the  Lefschetz  Center  for  Dynamical  Systems  for  its  invitation  to 
give  these  lectures,  and  to  the  audience  for  its  patience  and 
its  motivating  questions.  Special  thanks  are  due  Professor  H.  T. 
Banks  for  his  helpful  suggestions  and  comments  while  the  author 
was  preparing  these  notes. 
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CHAPTER  I 


Introduction  to  Immobilized  Enzyme  Systems 
I. 1.  What  are  enzymes? 

Enzymes  are  catalysts  of  biochemical  reactions 
In  the  simplest  case  typically,  one  reactant  S 
strate,  is  consumed  under  the  catalytic  action  of  an  enzyme  E, 
and  yields  a  product  P. 


n:  y 

ctions . 

,.  called  a  sufexV' 


(1.1) 


H 


-*>  P 


For  example,  in  such  a  reaction  one  might  have  E  =  glucose  oxidase, 
S  =  glucose,  P  =  gluconic  acid.  More  than  one  substrate  can  be 
involved  in  the  reaction.  A  second  substrate  is  often  called  a 
cosubstrate.  Similarly,  more  than  one  product  can  result  from 
the  reaction.  For  example,  with  S  =  uric  acid,  A  =  oxygen, 

E  *  uricase,  P  =  allantoin,  we  can  write 

E 

(1.2)  S  +  A  - —  *  P  +  other  products. 

Generally,  an  enzyme  catalyzes  only  a  single  type  of  reaction. 

Enzyme  kinetics  have  been  studied  since  at  least  Michaelis  and 
Menten  [5]  and  a  detailed  description  can  be  found  in  Dickson  and 
Webb  [6J. 


The  rate  expressions  for  enzyme  reactions  in  a  well-stirred  solution 
vary  greatly.  One  of  the  most  simple  is  the  monoenzyme  irrever- 
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sible  reaction  rate  expression 


(1.3) 


r 


[Si 

KS+[S] 


Here  Kg  -  Michaelis  constant,  characteristic  of  the  enzyme; 
VM  ■  maximal  reaction  rate,  proportional  to  enzyme 
concentration  in  the  solution;  and 
[S]  *  substrate  concentration. 


We  shall  deal  with  substrate  inhibited  phenomena  for  which  the 
rate  is  given  by 


(1.4) 


'M 


[S] 


K  ♦iS](l+J£I) 

5  ss 


Here  Kss 


inhibition  constant. 


In  the  case  of  the  reaction  (l .  2)  ,the  rate  expression  may  be 

taken  as 


(1.5) 


r  -  VM  - [S] - 

KA+[A1  K.+  tSKl+X^l) 
s  KSS 


Here  [A]  *  co  substrate  concentration  and 

■  Michaelis  constant  of  enzyme  E  for  cosubstrate  A. 


If  [A]  is  small  with  respect  to  [K.],  we  can  approximate  (1.5)  by 


In  the  following  we  shall  choose  Kg  as  the  unit  of  con 
centration,  and  use  the  dimensionless  quantities 


The  rates  of  reaction  (1 . 3) , (1. 4) , (1 . 5) ,  and  (1.6)  can  then  be  written 


For  a  detailed  discussion  of  enzyme  kinetics,  see  Banks  [7]  and 
Murray  [3]. 
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1.2.  What  are  immobilized  enzymes? 

In  living  cells,  enzyme  molecules  are  not  in  an  homogeneous 
solution,  but  are  linked  to  other  protein  molecules  within  a 
medium  that  is  heterogeneous  requiring  the  diffusion  of  metabolites 
(the  substrates  and  the  products).  In  other  words  there  exist 
gradients  of  concentration,  and  the  interaction  of  diffusion  and 
enzyme  reaction  must  be  taken  into  account  in  the  understanding 
of  enzyme  kinetics  in  living  cells.  Of  course,  the  phenomena  here 
are  much  more  complicated.  Other  phenomena  include  convection, 
the  possibility  of  electric  fields  due  to  ionized  species,  and 
the  co-occurence  of  many  enzyme  reactions,  all  of  which  make  the 
study  of  a  real  cell  very  complicated. 

Thus,  the  technology  of  artificial  enzyme  membranes  has 
been  developed  (Thomas  [8]):  this  enables  one  to  study,  in  a 
well-defined  context,  basic  phenomena  due  to  the  interaction  of 
diffusion  and  enzyme  reaction.  Additional  motivation  lies  in  the 
potential  use  of  artificial  enzyme  membranes  for  industrial, 
analytical  and  medical  applications.  We  refer  for  a  discussion 
of  this  aspect  to  Thomas  [9] . 

By  immobilized  enzyme  system  we  mean  artificial  membranes 
produced  by  linking  enzyme  molecules  to  inactive  protein  molecules, 
the  whole  forming  a  matrix  where  the  enzyme  is  uniformly  dis¬ 
tributed.  Such  a  membrane  is  a  translucent  film  of  about  50y 
thickness  which  is  easily  handled  in  the  laboratory. 


1.3.  The  "model  case'*. 

In  the  simplest  case  an  artificial  membrane  M  (slab 
geometry.  Fig. 1 . 1) separates  2  compartments  CQ  and  ,  con¬ 
taining  substrate  S  at  fixed  concentrations  Sq  and  s^.  S 
diffuses  into  M  and,  under  the  catalytic  action  of  E,  reacts: 
one  molecule  of  E  and  one  molecule  of  S  combine  to  form  a  molecule  of 
enzyme  -  substrate  ES,  which  itself  yields  one  molecule  of  product  P  and 
releases  one  molecule  of  enzyme  E:  schematically, 

(1  .12)  E  +  S?=*ES— s*>P  +  E. 


Thus,  the  net  action  of  the  enzyme  is  such  that  S  dis¬ 
appears  in  the  membrane  with  a  rate  given  by  (1.3). The  local  varia¬ 
tion  of  substrate  concentration  within  the  membrane  is  due  to 
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diffusion  and  reaction;  this  is  expressed  in  the  equation  of 
continuity : 


(1.13) 


0. 


Here  [S] 
t 


x 


.  3 

concentration  (moles  cm  ) 
time  (hour) 

2  - 1 

diffusion  coefficient  (cm  h  ) 

abscissa  in  the  direction  transverse  to  the  membrane (cm) . 


By  taking  as  new  units  Kg  for  concentrations,  L(thickness 

2 

of  the  membrane)  for  length  and  0(=L  /Dg) (time  lag)  for  times,  we 
obtain  the  equation 


(1.14) 


as 

at 


o 


0  <  x  <  1 ,  t>0 


where  the  single  parameter  a  is  of  importance  for  the  behavior  of 
the  membrane: 


(1.15)  a  =  jj-  (square  of  the  Thiele  modulus). 

S  S 

Equation(l . 14)must  be  supplemented  with  boundary  conditions 


s (0 , t)  =  s  0  ,  s ( 1 ,  t )  =  sx 


(1.16) 


and  initial  conditions;  for  example. 


(1.17)  s (x,0 )  =  0 

which  corresponds  to  a  membrane  initially  empty  of  any  substrate. 

The  proof  of  existence  and  uniqueness  of  a  (positive) 
solution  for  eqs .  (1 . 14)  ,  (1 . 16)  ,  (1 . 17)  ,  as  well  as  for  the  other 
equations  of  evolution  presented  in  this  chapter,  can  be  found  in 
Kernevez  and  Thomas  [10], 

Numerical  simulations  (using  finite  differences)  produce 
an  evolution  of  the  profile  of  concentration  tending  to  a  stable 
steady  state  (Fig. 1.2). 


Figure  1.2 
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The  main  interest  of  this  example  is  to  demonstrate  possibility 
of  a  fact  which  is  not  always  taken  into  account:  the  existence 
of  gradients  of  concentration  in  the  interior  of  the  membrane. 
Artificial  membranes  coupled  with  electrodes  (an  important  con¬ 
figuration  in  many  applications)  also  fall  within  this  example: 
a  membrane  is  attached  to  an  electrode  which  is  sensitive  to  the 
product  P  of  the  enzyme  reaction  (Fig.  1.3). 
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Figure  1.3 


Due  to  the  presence  of  the  electrode,  there  is  a  zero-flux  boundary  condition 
at  x  =  1/2  so  that  the  S  profiles  evolve  like  the  left  half  of  the  profiles 
in  Fig.  1.2,  whereas  the  P-profiles,  governed  by 


mi 


(1.18) 


(1.19) 


Cl . 20) 


D  .2 


p(0,t)  =  0  ,  -  0  for  x  =  j- 


P (x , 0)  =  0, 


evolve,  as  depicted  in  Fig. 1 . 3, towards  a  stable  steady  state, 
Of  course  in  this  example (1 . 16) has  to  be  replaced  by 


(1.21) 


s(0,t)  =  s. 


•  §7(7^  ■  ° 


From  the  measurement  of  the  electrode  one  can  get  the 
concentration  of  P  against  the  electrode,  and  hence  the  con¬ 
centration  of  S  at  x  =  0,  that  is,  in  the  solution  where  the 
electrode  is  immersed.  Due  to  the  specificity  of  the  enzyme,  the  electrode  is 
specific  for  one  particular  substrate  (glucose,  saccharose,  lactose, 
etc.).  It  is  thus  possible  to  monitor  continuously  the  concentra¬ 
tion  of  these  substances  in  media  like  the  blood  or  the  broth  of 
fermentations . 


1.4.  The  glucose  pump.  (Thomas  et.al-  [8]) 

In  this  example  there  is  a  spatial  distribution  of  the 
enzymes  and  the  interaction  of  diffusion  and  reaction  creates  a 
transport  of  S  against  a  concentration  gradient. 


The  membrane  (Fig. 1 . 4)consists  of  2  layers,  one  with  enzyme  E 
(hexokinase)  and  the  other  with  enzyme  E2  (phosphatase). 
transforms  S  (glucose)  into  P (glucose-6-phosphate)  and  E2 
transforms  P  into  S: 


(1.22) 


impermeable  barrier 
for  p 


The  membrane  separates  2  compartments  CQ  and  Cj,  with  sub¬ 
strate  concentrations  sQ  and  s1  (Sj^  _>  sQ) .  However,  the 
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system  pumps  glucose  from  CQ  to  C^.  This  is  seen  from  the 
form  of  the  S  concentration  profile  (Fig. 1.5). 


1 

I 

I 

I 

1 

I 


s 


S  is  consumed  for  0  <  x  <  j;  its  profile  is  convex.  S 
is  produced  in  the  region  j  <  x  <  1;  its  profile  is  concave. 

The  flux  of  S  at  x  *  0  (resp.  x  =  1)  being  given  by 
fresp.  -  It'1  ,t)),  we  see  that  S  is  entering  into  the  membrane 
at  x  =  0,  and  going  out  of  it  as  x  =  1.  The  result  is  an  active 
transport  of  glucose  from  CQ  to  C^;  this  has  been  observed 
experimentally,  and  also  is  found  when  one  solves  numerically  the 
corresponding  equations: 
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3t 


3x 


9s  ♦  r  .  o,  II  -  if$  .  r  -  0 


ax* 


r  "  0  X+s+p  for  0  <  x  <  j  and  r  -  -o  i|p  for  \ 


<  X  <  1 


(1.23W 


s(0,t)  *  s0,  s(l,t)  •=  Sj_,  g(0#t)  *  ||(l,t)  -  0 


s(x,0)  =  0  ,  p(x,0)  =  0, 


Indeed,  numerical  simulations  confirm  the  possibility  of  active 
transport  against  the  concentration  gradient  even  for  Sq  ■  1 
and  Sj  *  100. 


1.5.  Substrate  inhibited  kinetics. 

All  the  systems  to  be  studied  henceforth  will  be  endowed  with 
substrate  inhibited  kinetics: 

(i)  systems  for  which  multiple  steady  states  may  exist. 


f  ds 
3T  " 


a 


s 

1+s+ks  2 


0,  0  <  x  <  1,  t  >  0 


(1*24) 


< 


s (0 ,t) 


s0 


S(l,t)  »  Sj 


(  s(x,0)  ■  fixed. 


These  equations  represent  the  evolution  of  substrate  concentrations 
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in  a  membrane  separating  2  compartments  with  fixed  concentrations 
sQ  and  Sj  (Fig. 1.1).  We  have  not  specified  the  given  initial 
profile  of  concentration  x  -*■  s(x,0)  but,  as  we  shall  see  in 
chapter  2,  the  steady  state  will  depend  upon  it.  a  and  k 
have  been  defined  in  (1.15)  and  (1.9). 

(ii)  systems  for  which  periodic  oscillations  in  time 


may  exist. 


+  aa 


1+s+ks' 


0  <  x  <  1,  t  >  0 


(1.25) 


32a 

a — i  +  aa 


l+s+ks' 


s(0,t)  =  sQ,  s(l,t)  =  sx,  a(0,t)  =  aQ,  a(l,t)  =  a1, 


Here  the  rate  of  reaction  is  that  given  in  (1.11)  and 


(1.26) 


VM  L2  „  DA 

0  *1^  '  “ 


where  DA  is  the  diffusion  coefficient  for  A.  (Though  a  is 
the  concentration  of  a  cosubstrate  A,  we  denote  it  by  a  be¬ 
cause  A  acts  like  an  activator  in  the  reaction.  ) 

(iii)  the  S  -  A  system. 

We  consider  a  reaction  such  as  (1.2),  with  rate  (1.6) 
taking  place  in  an  active  layer  (i.e.  the  enzyme  is  immobilized 
within  it).  This  active  layer  is  thin  enough  to  be  considered  as 


L...  .. 
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a  planar  domain  ft  (Fig. 1.6)  where  the  concentrations  of  S 
and  A  depend  on  the  position  in  ft  but  not  on  the  position 
in  the  transverse  direction.  An  inactive  (i.e.  without  enzyme) 


reservoir 


inactive  layer 


active  layer 


i  i 


boundary  layer  separates  ft  from  an  outside  well-stirred 
reservoir  where  S  and  A  are  at  fixed  concentrations  SQ  and 
A q.  In  ft  3  phenomena  take  place:  transport  from  the  reservoir 
through  the  inactive  layer,  diffusion  and  reaction.  The  balance 
of  masses  for  S  can  be  written 
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Here  (resp.  L2)  is  the  thickness  of  the  active  (resp. 

inactive)  layer,  and  Dg  (resp.  Dg)  is  the  diffusion  coefficient 

of  S  in  the  active  (resp.  inactive)  layer. 

Taking  as  new  units  L  (diameter  of  D)  for  length, 

2 

9  (■  L  /Dg)  for  time  and  Kg  for  concentrations  we  obtain 


(1.28) 


”  -  As  +  aa  - - t- 

3t  l+s+ksZ 


-  Y(sn-s)  =  0. 


Here  A  is  the  2-dimensional  Laplacian,  0=^—0  and  y  -  0/0' 

i 

(0'  »  L1L2/Dg) .  Similarly  for  A  we  find 


(1.29) 


-  8Aa  +  oa 


l+s+ks' 


ay(a0-a)  =  0, 


D.  D. 

where  a  m  -  and  0  ■  are  the  ratios  of  the  diffusion  co- 

DS  ^ 

efficients  for  A  and  S  in  the  inactive  and  active  layers. 

Equations  (1.28)and(l.29)are  supplemented  with  zero  flux 
boundary  conditions: 


(l . 30) 


3s  _  />  j  3 a 
•k —  0  and  — 

3n  3n 


(|^  ■  n  •  7s,  n  being  the  unit  normal  to  3D  pointing  outwards). 
As  usual,  initial  conditions  for  s  and  a  have  to  be  specified. 
We  shall  deal  with  the  S  -  A  system  in  chapter  3  for  wave  front 
propagation,  and  in  chapter  4  for  pattern  formation.  There  we 
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shall  be  mainly  interested  in  its  steady  states.  The  S  -  A 
system  more  generally  can  be  the  modelling  of  diffusion  and  enzyme 
reaction  on  a  surface  separated  from  an  external  medium  by  a 
boundary  layer.  The  Laplace -Beltrami  operator  Ag  is  then  appro¬ 
priate  to  represent  diffusion  on  the  surface,  and  the  S  -  A 
system  equations  are 

If  -  V  +  aF(s’a)  *  Y(s0-s)  =  0 

(1.31)  <  -  BA^a  +  aF(s,a)  -  ay(aQ-a)  =  0 

,  F(s,a)  -  as/ (l+s+ks^) 

with  no  boundary  conditions  if  the  surface  is  closed. 

the 

Let  us  recall  the  definition  of/Laplace  Beltrami  operator: 
locally  the  surface  (think  of  a  cucumber  surface.')  can  be  mapped 
into  the  interior  of  a  square:  -1  <  u^<  +1,  -1  <  ^  <  +1  (Fig. 1.7) 


Define 


Then 


(-AB<|>)i|ids 


As  we  shall  see  in  chapter  4,  it  is  no  more  difficult  to 
solve  numerically  the  S  -  A  system  on  a  cucumber  surface  than 
on  a  planar  surface  by  using  the  finite  element  method.  However, 
there  we  shall  have  in  mind  the  surface  of  an  egg  rather  than  the 
surface  of  a  cucumber! 
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CHAPTER  II 

Multiple  Steady  States  and  Hysteresis 
2.1.  The  0 -dimensional  case. 

The  equation 


(2.1)  Jt  3  si?  "  s  '  PF(S) 

represents  the  rate  of  change  of  the  substrate  concentration  s 
in  a  well-stirred  cell  (Fig. 2.1)  where  a  substrate  inhibited  reaction  takes 
place : 

(2.2)  F(s)  *  - 

l+s+ksz 


so 

reservoir 

-X- 

cell 

membrane 


Figure  2.1 


The  cell  is  separated  from  a  well  stirred  reservoir  with  fixed 
concentration  sQ  by  an  inactive  (i.e.  without  enzyme)  membrane 
transporting  the  substrate  by  diffusion  from  the  reservoir  to  the 


T 

L 
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cell  with  a  flux  sQ  -  s.  In  this  system  diffusion  and  (enzyme) 
reaction  take  place  in  distinct  spatial  regions;  however,  it  is 
one  of  the  simplest  imaginable  systems  in  which  the  coupling  of 
reaction  and  diffusion  yields  a  phenomenon  involving  multiple 
steady  states;  in  this  case  governed  by  solutions  of 


(2.3) 


sn  -  s  =  PF(s). 


In  fact,  it  is  also  an  open  system,  due  to  the  transport  into 
the  system  with  a  flux  sQ  -  s.  Without  this  transport  from 
outside,  the  system  would  evolve,  according  to  the  kinetic 
equation 


(2.4) 


|f  =  -  PF(s), 


to  the  equilibrium  state  s  =  0.  But  equation  (2.2)  possesses 
solutions  "far  from  this  equilibrium",  simply  obtained  by  cutting 
the  curve  y  ■  pF(s)  by  the  straight  line  y  ■  sQ  -  s  (Fig. 2. 2), 


y  -  pF (s) 


S1  s2 


Figure  2.2  -  solutions  of  (2.2). 


If  we  fix  p  in  (2.3)  and  wish  to  get  s  as  a  (possibly 

multivalued)  function  of  Sq,  we  can  consider  the  graphs  of 

s  s0  =  s  +  pF(s)  (Fig.  2.3)  and  of  s  -*■  F'(s)  where  it  is 

easily  seen  that  if  p  >  -1/min  F'(s),  there  are  2  values  s* 

ds0 

and  s**  of  s  such  that  *  1  +  pF'(s)  =  0,  and 

Sq  *  s  +  pF(s)  admits  relative  extrema. 


Figure  2.3 

In  that  case  we  shall  have  the  representations  depicted  in  Fig. 
2.4  for  s  as  a  function  of  s0  and  for  the  reaction  rate  of 
the  system: 


*  pF(s)  . 


So  -  S 


21 


All  of  the  above  is  easily  seen,  as  is  the  stability  of  the 
state  s  with  respect  to  the  dynamical  system  (2.1)  which  changes 
when  1  ♦  pF'(s)  changes  sign,  i.e.  at  the  turning  points 
s  ■  s*  and  s  ■  s**.  However,  the  behavior  of  the  simple  system 
shown  in  Fig.  2.4  will  extend  to  more  complex  cases,  e.g.  dis¬ 
tributed  systems. 

The  O-dimensional  case:  2  cells.  (Bunow  et. al.  [11]) 

We  consider  now  2  cells,  (see  Fig<  2.5),  where  the  substrate 
concentrations  Sj  and  S£  are  governed  by 


. . -  ■ . . — ~ 
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(2.6) 


ar* +  2si  '  s2  +  pF(si)  -  s( 

ds2 

at-  '  S1  *  2s2  +  pF(s2)  *  s. 


Figure  2.5 


^  branch  of  solutions  of  the  steady  state  equati 


(2.7) 


2sl  '  s2  +  pp(si)  “  sQ 


~S1  +  2s2  +  pp(s2)  *  s0 


consists  of  "symmetric  solutions"  sx  -  s2  =  s,  with  (2.3).  But 
there  may  exist  "non  symmetric"  steady  states  (sx  j  s2),  as 
follows  from  the  stability  study  of  the  symmetric  states 
(S1 » s2^  *  (s,s):  the  linearized  operator 


' 


2+pF ' (Sj) 


-1 


2+pF'Cs) 


-1 


t 

- 

1  L  -1  2  +  pF ' (s?) . 

-1 

2+pF' (s) 

has  eigenvalues  »  1  +  pF'(s)  and  X2  *  3  +  pF'(s).  Thus  in 
addition  to  the  turning  points  (sg,s*,s*)  and  (sg*,s**,^*) 
corresponding  to  X^  =  0,  one  may  have,  if  p  is  large  enough, 
X2  =  0  (Fig.  2.6),  and  consequently  a  bifurcation  point.  As  the 
reaction  rate  is  the  same  for 2  "mirror  states"  (s1,s2)  and  (s2,s1) 


1 


the  reaction  rate  of  the  twobifurcated  states  will  be  represented  by 
only  one  curve  (Fig.  2.7). 
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Figure  2.7 


One  interesting  feature  of  Fig.  2.7  is  that  stable  non-uniform 
steady  states  (sj  f  s2)  exist  on  these  bifurcated  branches. 

A  nontrivial  problem  is  how  to  attain  such  a  N.U.S.S. 
(Sq.s^s^  for  the  dynamical  system 


or  for  a  sequence  of  steady  states  (s^,s2) 


(2.9) 


|  2s1  -  s 2  +  pFfSj)  =  v1(t) 
-s1  +  2s2  +  P F (s 2 )  -  v2(t) 


In  both  cases  we  should  have,  at  some  given  time  T, 


(2.10)  s1(T)  =  s1,s2(T)  =  s2,v1(T)  =  v2(T)  =  s( 


Imperfection 

If,  instead  of  (2.7),  one  has 


(2.10a) 


2s1  -  s2  +  pl' (s -| )  =  s 


1J  J0 


-s1  +  2s2  +  pF(s,)  =  sn  +  e 


2J  * 0 


then  we  have  no  more  a  bifurcation,  but,  in  Fig.  2.7,  a  "lemniscate 
detaches  from  the  main  branch  (Fig.  2.8). 


t 

e  =  0 


s'-s 

X  /  / 

*  ^ 


/  e  f  0 


Fieure  2.8 
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This  lemniscate  comes  in  part  from  one  of  the  bifurcated  branches 
and  in  part  from  the  arc  of  symmetric  states.  It  is  an  example 
of "imperfection"  (Fig.  2.9). 


critical  value 

when  *  ’<=  1^4 
e  =  0—  — <  — 

N\ 

\ 

5  2* 


Figure  2.9 


2.2.  1-dimensional  case. 

We  consider  next  the  system  governed  by  equations  (1.24) 
of  chapter  1,  and  we  wish  to  study  the  steady  states.  They  satisfy 
the  2 -point  boundary  value  problem 

,2C 

4-4  +  oF(s)  *  0  ,  0  <  x  <  1, 

dxZ 

F(s)  =  s/(l+s+ks2) 
s (0)  *  s(l)  *  s0. 


(2.11) 


Here  we  take  s(0)  =  s(l)  in  order  to  simplify  the  analysis. 
Problem  (2.11)  has  a  first  integral 


(2.12) 


(x)  +  2<jG  (s  (x)  ) 


-fcG(u) 


where  G  is  the  indefinite  integral  of  F(G'=F)  and  p  =  s(l/2). 

(Since  s"  >  0,  s  is  a  convex  function  of  x,  and  since 
s(0)  =  s(l)  the  s  profile  is  symmetric  with  respect  to  x  *  j.) 
Moreover,  if  we  restrict  our  considerations  to  ^  <  x  <  1, 
s  '  (x)  >0  so  that 


=  (2a)1/2(G(s)-G(y))1/2 


and 


(2.13) 


dS 

(G  (C)“G(vi)  ) 1/2 


The  solution  of  (2.11)  is  equivalent  to  the  solution  of  (2.13) 
together  with  the  condition  requiring  s(l)  =  sQ: 


(2.14) 


,0,1/2  m  [S°  dg 

7  \  (G(5)-G(y))1/2 


In  fact  (2.11)  has  as  many  solutions  as  we  can  find  several  values 


of  M  satisfying  (2.14),  or 


where 


It  can  be  shown  (Kernevez  [10])  that  for  Sq  large  enough 
the  graph  of  f  has  the  shape  indicated  in  Fig.  2.10,  so  that 
there  are  3  solutions  to  (2.15)  if  a  is  in  a  suitable  interval. 


The  upper  and  lower  solutions  can  be  easily  obtained  numerically 
by  using  the  following  scheme: 

(i)  start  from  s^(x)  =  Sq  for  the  upper  solution 
and  s®(x)  *  0  for  the  lower  solution: 

(ii)  define  s^+1^(x)  from  s  ^  (x)  by 


and  iterate  until 


max  |  s 


One  thus  obtains  two  stable  steady  states  for  the 


same  boundary  value  s 


Figure  2.11 
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Fot  a  detailed  study  of  the  multiple  steady  states  in  the  one¬ 
dimensional  case  we  refer  to  Brauner  and  Nicolaenko  [12].  The 
curve  which  represents  yC-sfj))  as  a  function  of  s0(*s(0)=s(l)) 
is  similar  to  the  one  in  Fig.  4.  Hysteresis  may  be  observed 
for  a  given  boundary  value  Sg  the  system  can  be  in  one  of 
two  possible  stable  states  depending  on  its  past  history.  Thus 
we  have  a  simple  model  of  a  system  with  memory  (see  Fig.  2.12). 


Figure  2,12 


2.3.  2-dimensional  case. 


We  are  interested  in  the  steady  states  defined  by 
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I 


i 

I 


(2.19) 


-As  +  oF(s)  =  0 
F(s)  =  s/Cl+s+ks2) 


in  ft 


s0 


2  3 

where  ft  is  a  bounded  set  in  R  or  R  .  Brauner  and  Nicolaenko 
[13]  have  demonstrated  the  possibility  of  multiple  solutions  for 

(2.19) .  Their  proof  uses  a  regular  perturbation  method:  equation 

(2.19)  can  be  rewritten 


(2.20) 


-Au  + 


u 


7 - 7 — 7 

£  ^+8eu+kB^uz 


s 


1 


0 


2  1 

by  defining  u  =•  s/s0,  e  *  -  ,  $  =  sQe. 

We  have  already  made  the  remark  that  in  order  to  obtain 

multiple  steady  states  we  must  have  o  and  Sq  sufficiently 

1  B 

large.  Therefore,  the  hypothesis  a  *  — j  and  so  “  f"  seems 
reasonable,  and  in  fact  in  the  one-dimensional  case,  typical  values 
for  which  there  are  multiple  states  are  o  -  1200,  s0  =  72, 
corresponding  to  e~ 57*  6  “  2.2. 

Brauner  and  Nicolaenko  prove  the  following: 

(i)  there  exists  a  branch  of  maximal  solutions  to 
f  -Au  +  — i-  =  0 

J  kBzu 


1 


with  a  turning  point  T  (Fig.  2.13), 
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(ii)  for  e  small  enough,  the  same  is  true  for  (2.20), 
with  turning  point  T^, 

(iii)  there  is  at  least  one  other  turning  point  T2 
for  (2.20). 

Because  of  the  difficulty  in  obtaining  theoretical  results 
on  the  qualitative  behavior  of  the  solution  of  equations  such  as 
(2.19),  it  is  interesting  to  develop  efficient  algorithms  to 
follow  the  branch  of  solutions  when  a  parameter  like  s^  varies. 
Such  an  algorithm  is  described  In  the  next  section. 


Figure  2.13 


Graph  of  u(A),  A  being  any  point  in  n. 


2.4.  Numerical  analysis 


A  finite  element  discretization  of  (2.19)  gives  an 
equation  of  the  form 


is  the  vector  of  (unknown)  nodal 


Here  X  =  s 


so  that  in  fact  (2.21)  is  a  set  of  N  (non-linear)  equations  of 
N  unknowns  depending  upon  the  parameter  X. 

The  method  of  Kubicek  [14]  consists  of  solving  the  con¬ 


strained  system 


where  t  is  the  arclength  on  the  curve  of  solutions.  This  is 
equivalent  to  solving  the  differential  equations 


34 


for  which  we  can  employ  any  numerical  method  (Euler,  Adams- 

dX 

3t  • 


Bashforth,  ...)  if  we  are  able  to  calculate  ^  and  ^ 


We  can  proceed  in  the  following  way: 

(i)  Solve 


(2.2S) 


3F 

3T 


(ii)  We  know  that  gx  =  X  •  Replace  in  the  2nd 


equation  of  (2.24): 
2 


(2.26)  (|y|£+l)(§£)  =  1, 


dX 

3t 


±1 


/pi 


y|*i 


Suppose, for  example,  that  we  are  on  an  arc  of  the  curve 

X 

dt 


where  44  >  0 .  Then  'iQr  is  known  and  hence 


dX 

3T 


dx 

at  • 


(iii)  The  application  of  Euler  method  for  instance  will 
make  the  point  (x,X)  move  from  to  M2  (Fig.  2.14). 


Figure  2.14 
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(iv)  The  Newton-Raphson  method  or  the  secant  method 
applied  to  (2.21)  will  then  give  on  the  curve.  The  drawback 

of  this  method  is  that  when  one  approaches  a  turning  point,  the 
jacobian  -g^,  which  can  be  monitored  by  its  determinant,  has  a 
tendancy  to  become  singular.  As  we  need  to  solve  (2.25)  and 
linear  systems  with  the  same  matrix  for  Newton's  method,  we 

have  to  check  whether  its  determinant  is  not  too  small.  If  it 
is,  then  we  have  to  interchange  a  spatial  variable  (xN,  for  instance) 
and  X:  the  first  equation  of  (2.24)  can  be  written 


(2.27) 


Here  x'  * 


(2.28) 


3F ' 
3x ' 

3FN 


3F 

3x 


3F 


1 

'dx'  1 

'  3F'  " 

ar~ 

+ 

3X 

dxN 

3FN 

L~dt . 

.  3X  - 

dX 

ar 


txlx2  * 

•  *XN- 

1]T,  F* 

=  [Fr. 

•fn-i 

]T  • 

3F ' 

dx' 

+  3FJ_ 

dxN 

3F ' 

dX  _ 

W1- 

eft- 

3xn 

“at 

3  X 

at  " 

3Fn 

dx' 

,  3Fn 

dxN 

3Fn 

dx  . 

a’x' 

at- 

“at 

3X 

at  " 

or 


(2.29) 


3F' 

3xt 


3Fn 


3F ' 
3  X 


3Fn 

nr 


r  dx  *  i 

‘  3F '  ' 

+ 

3xn 

dx 

sfn 

J 

.eft  . 

J 

dxN 

w 
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Solving  (2.29)  together  with  the  2nd  equation  of  (2.24), 
or  solving  (2.24),  are  similar  problems.  We  essentially  have  in 
each  case  to  solve  a  system  of  the  form 


A11 

A12 

"yl" 

'fl* 

.  A21 

\ 

CM 

N 

< 

_y2 . 

.f2. 

where  An  is  always  the  same  (N-l)  x  (N-l)  matrix  ,  A.  7 

is  (N-l)  x  1,  A22  is  lx  (N-l)  and  A22  €  IR.  We  want  to  take 
advantage  of  the  fact  that  is  symmetric  and  sparse  and  use 

skyline  storage  for  it  together  with  an  efficient  solver  (Bathe 
and  Wilson  [15]).  Suppose  that  we  have  already  made  the  triangu- 
larization  of  An  :  An  =  LDLT.  We  observe  that 

Allyl  +  A12y2  fl  *  yl  Allfl  '  AllA12y 2  * 

But  Cj  *  Ai]^i  and  c2  *  AllAl2  can  be  foun<*  by  solving 

(2.31)  AnCl  - 

and 


(2.32) 


A11C2  "  A12  * 


! 
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Now  replace  by  Cj  -  Y2C2  *n  t*ie  2nd  equation  of 

(2.30)  ( ^2 2  *  ^2  an<*  ^2  are  numbers)  : 


A21c1  '  y2A21c2  +  A22y2  "  f2 


Hence, 


(2.33) 


2.5.  Optimum  design. 

We  shall  investigate  the  dependence  of  solutions  of  (2.19) 
upon  the  geometrical  domain  q,  using  methods  of  Murat  and  Simon 
[16],  Mignot,  Murat  and  Puel  [17], 

For  a  given  domain  0  and  a  given  value  of  a  we  know 
that  the  activity 


(2.34) 


a 


(  i 1 
J  3  n 


dt 


F ( s ) dx 


varies  with  sQ  as  indicated  in  figure  2. IS.  We  are  interested 
in  the  value  of  sQ  for  which  the  activity  is  maximum. 


1 

I 


In  the  following  sj  and  s*  will  denote  this  sQ  (a 
number)  and  the  corresponding  state 

(2.35) 

(2.36) 

(2.37)  -AS  +  oF’(s*)S  *  0,  S|  *  1. 

Ir 

Now  we  keep  a  fixed  but  allow  ft  to  vary  in 


-As*  +  aF(s*)  =  0,  s' 


%  *  0  or  {  H  dr  =  0  "here 
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I 

1 

I 

I 

1 


(2.38)  =  {ft|ft  *  T (Q)  ,  T  €&2  } 

? 

where  Q  is  a  given  domain  with  C  boundary  and 


(2.39) 


%2  *  {T  | T- 1  and  T'1-!  G  (Rn  ,IRn)  } ,  n  =  2  or  3. 


Moreover,  we  impose  the  condition 


(2.40) 


meas(Q)  =  k. 


Since  the  total  quantity  of  enzyme  in  the  membrane  is  proportional 
to  o  meas(ft),  this  quantity  is  constant  when  the  shape  of  ft  varies. 
We  are  interested  in  observing  the  influence  of  the  shape  of  ft  on 
the  maximum  activity  a*.  This  is  not  motivated  by  any  direct  in¬ 
dustrial  requirement,  but  is  motivated  by  our  desire  to  study  the 
interaction  of  diffusion  and  reaction  in  simple  artificial  enzyme 


systems . 

We  shall  need 

the  2  formulas 

(2.41) 

tff(f  fftdx)T  " 

f  3£n 

J  SIT  T  dx  *  J 

» 

fn  n • xdr 

ft  7 

ft 

r 

(2.42) 

Sn(J  gftdr)T  " 

j 

[  (stt  ♦ 

r  r  r 


Here  ffi  and  gfl  are  functions  depending  on  the  parameter  ft 
and  H  is  the  mean  curvature  of  the  variety  3ft  ■  r. 


(me as  ft) 


n*T  dr.  If  meas  (ft)  *  k 


where 


Therefore 


F ' (s  *) S*dx 


S*  satisfies  the  conditions  obtained  by  differentiating  equations 
(2.35)  : 


The  boundary  condition  (2.45)  can  be  obtained  by  the  following 
trick.  Let  a  be  any  function  in  C°°(Rn).  Then 


(2.46) 


s*  *  sn  on  r 


|  (s*-sj) 


adr  =  0  Va  G  C*(IRn) 


Apply  (2.42) 


!ii(f  (s-s*)adr)T).  |  (s*-sj)  dr 
r  r 

*  f  [Hra  *  1*  *  H(s*-sJ)]  n-Tdr 


s*  -  s0  =  0  on  r 


0  ■  J  (S*-SQ)adr  +  j  a  n’xdr 

r  r 


whence  (2.45).  Now  let  us,  at  last,  consider  (2.36),  (2.37),  (2.43), 

* 

(2.44)  and  (2.45),  and  introduce  u  =  S*  -  SnS: 


a  F'(s*)S*dx 


3S*  Ar  q*  f  3S 

7n“  dr  '  S0  J  ITi  dr 


(2.47) 


and  u  satisfies 


(2.48) 


•Au  ♦  oF’ (s*)u  *  0  u 


In  conclusion  for  a  given  ft  the  maximum  activity 
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is  obtained  by  solving  (2.35),  (2.36),  (2.37),  and  the  deformation 
of  ft  in  order  to  increase  this  maximum  activity  (while  meas(Q) 
and  consequently  the  total  quantity  of  enzyme  remain  fixed)  is 
given  by  (2.47)  and  (2.48). 
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CHAPTER  III 

Oscillations  and  Wave  Front  Propagation 
3.1.  Oscillations:  O-dimensional  case 

One  of  the  simplest  situations  in  which  oscillations 
arise  is  when  a  stable  steady  state  looses  its  stability  by  under 
going  a  Hopf  bifurcation.  The  steady  state  becomes  unstable  and 
gives  rise  to  a  branch  of  periodic  solutions  (Fig.  3.1). 


— — »  stable  steady  state 

—  —  —  —  —  unstable  steady  state 

stable  periodic  solution 
ooooooooeoooooooo  unstable  periodic  solution 

| | s | |  =  amplitude  of  s,  state  of  the  system 
X  *  bifurcation  parameter 

Figure  3.1 


I 
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In  particular  if  the  system  is  governed  by  a  pair  of  coupled 
ordinary  differential  equations 


(3.1) 


(  ds 
at 


f (s,a) 


da 

3t 


g(s,a)  , 


and  if  the  "particle"  (s(t),a(t))  remains  in  a  rectangle 
0  £  s  £  Sq,  ®£a£a0»  then  we  satisfy  the  conditions  of 
Poincar^-Bendixson  theorem:  particle  lying  in  a  rectangle,  existence 
of  an  unstable  equilibrium  point.  Therefore,  there  exists  a  stable 
limit  cycle.  This  occurs  in  the  S  -  A  system  (Ch.  1,  §5). 


(3.2) 


ds 

cfit  =  so  ’  s  '  paF(s) 

at  ”  °(ao'a)  "  PaF(s) 

F(s)  *  s/(l+s+ks2)  ,  k  >  0, 


for  well  chosen  values  of  the  parameters  sQ,aQ,a  and  p  (k  is 
given  by  the  enzyme  kinetics) . 

In  this  kind  of  problem  it  is  very  useful  to  consider  the 
phase  plane  (s,a)  (Fig.  3.2)  where  the  null-cline  curves  f  ■  0 
and  g  ■  0 
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tangle  [0,Sq]  *  [0,aQ]  (where  the  trajectory  remains)  in  regions 
where  j|-  and  are  either  positive  or  negative,  whence  the 

shape  of  the  limit  cycle  (Fig.  3.2).  Besides  this  very  simple  par¬ 
ticular  case  one  could  be  interested  in  finding  parameter  values 
for  which  a  given  system  might  oscillate.  For  instance,  the  model¬ 
ing  of  particles  of  enzyme  membrane  in  a  CSTR  (continuously  stirred 
tank  reactor -see  Fig.3.3)can  be  done  by  assuming  that  everything 
occurs  as  if  there  were  2  cells,  the  CSTR  and  the  interior  of  the 
particles,  communicating  by  diffusion  through  the  particles.  (Fig. 
3.4) 


s  ,a 

outlet 

small  particles 
of  enzyme  membrane 

Figure  3.3  Figure  3.4 


S0’a0 

s ,  a 

O  O 

->■ 

O 

inlet 

V 

[S0]D 


enzyme  reactcj 
(interior  yfr' 
the  particles) 


inactive 
membrane 
(diffusion  only) 


The  model  is 


vl^Sl  =  D(Sq- [S] ) 


D, 


„ ts]-[s1] 
— E - 


(3.4) 
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Here  V  =  CSTR  volume,  [S]  =  S  concentration  in  CSTR,  D  =  flow 
rate  in  CSTR,  L  (resp  £)  =  thickness  (resp.  area)  of  the  membrane, 
Upon  normilizing  variables  in  (3.4),  we  obtain 


(3.5) 


ds  1  ,  .  1  ,  , 

at  =  e(so"s)  '  rr(s_si) 


Here  s  =  [S]/Kg,  sQ  =  [Sg]/Kg,  e  =  V/D,e'  =  VL/(£Dg).  We  have, 
similarly,  for  A  in  the  CSTR: 


(3.6) 


at  '  tCa0'a)  '  tT’(a"al) 


u  » 

where  g  =  jj--,  where  is  the  diffusion  coefficient  in  the  membrane.  In  the 
S 

cell  corresponding  to  the  enzyme  membrane  particles  we  have 


[S]  -  [S,] 

l - £ -  "  vr(  [S]  ,  [A]  )  , 


d[Sl] 

V  -4-  -  Ds  l 


where  v  is  the  volume  and  r  the  rate  of  reaction,  or 


(3.7) 


ds  i  . 

3t~  =  7T  (s"sl)  ‘  F(si»ai) 

91 


where 


9i  =  xnr  and  F(si,ai)  =  h 


Similarly, 


(3.8) 


^al  1 

ar-  =  (a_ai)  ■  F (si » ai) 


1 


The  result  of  this  modelling  is  a  set  of  four  O.D.E.'s 
(3.5),  (3.6),  (3.7),  and  (3.8),  including  six  parameters 

t  i 

s0,a0,e,9’ ,61,6.  The  parameters  9,9'  and  9^  are  characteris¬ 
tic  times,  9  for  the  flow,  the  two  others  for  diffusion.  Now 
it  is  not  evident  that  given  the  rate  of  reaction  F(s,a)  by  a 
set  of  measurements  (Fig.  3.5),  we  can  choose  the  six  parameters 
so  that  a  stable  periodic  solution  arises  as  a  consequence  of 
Hopf  bifurcation.  We  shall  have  to  use  numerical  analysis. 


Figure  3.5 

? 

First,  we  identify  the  parameters  9,9'  and  9^^  from 
experiments  without  reaction  (that  is  with  particles  of  membrane 
but  without  immobilized  enzyme  within  them).  Then,  we  explore  all 
the  possible  steady  states  for  the  possibility  of  a  Hopf  bifurca¬ 
tion: 
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(i)  Choose  (e,s,,a,)  and  calculate  s  and  a  from 

ds,  da, 

(3.7)  and  (3.8)  with  =  0,  and  sQ  and 

aQ  from  (3.5)  and  (3.6)  with  =  0. 

(ii)  Calculate  the  linearized  operator  of  (3.5),  (3.6), 

(3.7) ,  (3.8) 


and  calculate  its  eigenvalues. 

If,  all  these  eigenvalues  being  in  the  left  half  plane, 
it  occurs  that  a  pair  of  complex  conjugate  eigenvalues  crosses 
the  imaginary  axis,  then  we  can  hope  that  a  periodic  solution 
will  arise. 

Now,  we  have  to  make  this  calculation  for  all  the  triples 
(a.s-pa^)  in  some  region  of  IR  .  This  is  time  consuming  but  it 
is  practically  impossible  to  experimentally  find  the  desired 


.  _  . 
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values  of  0,sn  and  an  for  which  the  system  oscillates. 


3.2.  Oscillations:  1-dimensional  case 

Equations  (3.2)  represent  the  interaction  of  diffusion 
(Sq-s  and  ag-a)  and  reaction  (paF(s))  in  a  system  where 
these  two  phenomena  take  place  in  separated  spatial  regions. 

In  the  case  of  a  distributed  system  such  as  an  immobilized  enzyme 
membrane, diffusion  and  reaction  are  governed  by  similar  equations 

but  in  this  case  the  diffusion  terms  sQ  -  s  and  a(ag-a)  are 

i  ,  ,  3Zs  ,  3Za 

replaced  by  — j  and  a — 7  : 


If  -  7T  +  oaF(s)  =  0 


0<x<l,t>0 


(3.9) 


+  aaF(s)  =  0 

3t  3xZ 


F (s)  =  s/  (1  +  s+ks^) 
s  (0  ,t)  =  s  (1 ,  t)  =  s 0 


a(0,t)  =  a(l,t)  =  af 


This  system  can  oscillate.  For  instance,  we  choose  Sq  and  ag 
such  that  Sq  -  aaQ  =  0.  Now  if  we  look  for  steady  state  solutions 
of  (3.9) 


(3.10) 


d‘s  + 

37 

jaF(s)  =  0 

s (0)  -  s(l) 

+ 

dxZ 

oaF(s)  =  0 

a(0)  -  a(l) 

/ 


This  latter  equation  admits  one  solution  s  since  s  -*■  - 

1+s+ks 

is  monotone.  The  corresponding  steady  state  for  (3.9), 
s  =  s(x),  a  *  a(x),  can  be  analyzed  for  stability  by  calculat 
ing  the  eigenvalues  of  the  linearized  operator 


The  pair  of  complex  conjugate  leading  eigenvalues  behaves,  when 
o  varies,  as  depicted  in  Figure  3.6. 


Plot  of  the  two  leading  eigenvalues  X  and  X  when  a  varies 


£ 


We  see  that  when  a  enters  into  the  interval  (o^c^),  the 
state  (s,a)  loses  its  stability  and  a  family  of  periodic 
solutions  arise  by  Hopf  bifurcation.  There  are  two  possibilities 
(Figure  3.7). 


subcritical  bifurcation  supercritical  bifurcation 

Figure  3.7 
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From  a  numerical  point  of  view,  if  we  solve  the  equations  (3.9) 
for  a  slightly  larger  than  Oj,  we  find  a  periodic  oscillation 
of  the  s  and  a  profiles  of  concentration,  but  we  cannot 
decide  whether  the  bifurcation  is  sub-  or  super-critical.  Further 
calculations  would  be  necessary  to  check  that. 


3.3.  Wave  front  propagation 

We  consider  here  the  S 


£ 


-  1-4  +  Y[F(s,a)  - 


(3.13)  < 


"  8— j  *  Y[F(s,a)  - 
3xz 

F(s,a)  *  pas/(l+s+ks2) 
3s  _  3a 


3x  3x 


0  at  x 


-  A  system  equations  (Ch.I,  §5) 

(s0-s)]  =0  0  <  x  <  1 ,  t  >  0 , 

<*(a0-a)]  =  0 

and  x  *  1. 


We  suppose  that,  in  the  (s,a)  plane,  the  null-clines 


(3.14) 


J  f(s,a)  -  F(s,a)  -  (sQ-s)  »  0 
.  g(s,a)  -  F(s,a)  -  o(aQ-a)  -  0 


are  as  indicated  in  Figure  3.8. 


g  >  0 


They  possess  only  one  common  point  (§,£),  which  is  situated 
near  the  top  of  the  curve  f  =  0,  and  which  is  stable  for  the 
dynamical  system 


(3.15) 


In  OTder  to  understand  the  behavior  of  the  distributed 
system  (3.13),  let  us  first  describe  the  trajectory  of  a  particle 
subject  to  (3.15)  when  it  starts  from  Pn(Fig.  3.9). 


Figure  3.9 


Due  to  the  sign  of  and  in  (3.15),  in  the  regions 

delineated  by  the  isoclines,  this  trajectory  is  as  shown  in 
Figure  3.9:  going  South-West  in  a  first  phase  (fast  phase), 
going  South-East  in  a  second  phase  (slow  phase) ,  North-East  in 
a  fast  third  phase  and  North-West  towards  the  equilibrium  point 
E  in  a  last  phase  which  is  slow.  The  fact  that  some  phases 
are  fast  and  others  slow  is  due  to  the  respective  order  of  magni¬ 
tude  of  f  and  g  in  these  phases  of  the  movement. 

Now  in  the  case  of  the  distributed  dynamical  system  (3.13), 
instead  of  one  particle  we  have  a  continuum  of  them,  each  one 
corresponding  to  a  value  of  x,  0  <  x  <  1.  In  addition  to  the 
dynamics  due  to  the  terms  f  and  g,  these  particles  are  influ- 


i'W''-  '  ~  ■  ■ 
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enced  by  their  neighbors  by  diffusion.  In  fact  the  state  of 
the  system  at  a  given  time  t  can  be  conveniently  represented 
by  a  "snake"  (Fig.  3.10)  which  is  the  line  (s(x) ,a(x))g<x<1.  The 
"head"  of  the  snake  is  at  x  -  0  and  its  "tail"  at  x  =  1. 


0  s  (x) 


Let  us  start  with  the  system  (3.13)  at  the  (diffusionally 
stable)  steady  state  (s,a).  The  "snake"  is  concentrated  into 
the  point  (S,a)  and  the  S  and  A  profiles  are  as  shown  in 
Figure  3.11  a)  and  b) .  Then,  during  some  time  interval  (0,t), 
we  suppose  that  Sq,  instead  of  having  a  constant  value  s0, 
is  perturbed  near  x  *  0  (Figure  3.11  c)).  The  consequence  is 
to  modify  the  profiles  of  S  and  A,  mainly  that  of  S,  together 
with  the  shape  of  the  snake  in  the  phase  plane  (Figure  3.12). 


After  the  S  profile  has  been  so  perturbed,  we  can  restore  sQ 
to  a  uniform  value  Sq  for  0  <  x  <  1,  the  front  of  S  con¬ 
centration  between  and  x2  is  going  to  propagate  rapidly 

towards  x  =  1  (Figure  3.13). 


58 


When  the  front  attains  x  =  1 ,  it  drops  down  to  a  small  value 
of  S  concentration  and  from  now  on  all  the  s  and  a  con¬ 
centration  profiles  will  be  flat  so  that,  after  this  (fast) 
first  phase,  the  representation  of  the  "snake"  will  be  only  one 
point  in  the  (s,a)  phase  plane,  and  this  point  will  evolve 
like  the  particle  on  the  trajectory  of  Figure  3.9. 


4tl1  phase 


If  we  consider  the  signal  as  it  is  perceived  at  the  extremity 
x  ■  1,  s  »  s(l,t) ,  we  have  the  one  represented  in  Figure  3.15: 


_ 4th  phase: 

/  x:  slow  recovery 

L 

Lr — end  of  1  phase: 

L-  3r  phase: 

1  fast  drop 

f~  fast  burst 

\  J 

with  overshoot 

l+_a - - - z_ 

t 

\ - duration  2nd  phase 

of  the  perturbation 


Figure  3.15 
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CHAPTER  IV 

Pattern  Formation 
4.1.  Introduction 

Morphogenesis  in  insects  is  one  of  the  best-known  phenomena 
of  pattern  formation  in  embryos.  This  is  due  to  the  fact  that 
Drosophila  flies, for  example,  have  been  extensively  studied  for 
several  decades  by  many  zoologists.  They  have  observed  the  follow¬ 
ing:  parts  of  the  adult  epidermis,  such  as  a  wing,  a  leg,  a 
haltere,  the  eye-antenna  part  or  the  genital  part,  are  formed  of 
cells  which  are  differentiated  but  at  the  beginning  of  the  develop¬ 
ment  were  undifferentiated  within  what  is  called  an  imaginal  disk. 
In  fact  in  early  stages  of  development  an  imaginal  disk  is  formed 
of  very  few  cells,  say  40.  Marker  experiments  show  that  even  if 
at  the  beginning  all  the  cells  of  a  given  imaginal  disk  have  the 
same  potentialities,  after  some  time  there  is  a  difference  between 
their  commitments.  Take  the  wing  disc, for  example.  Cells  first 
differentiate  between  those  which  will  give  the  anterior  part  of 
the  wing,  and  those  which  will  give  the  posterior  part.  This  can 
be  seen  by  marking  a  cell  after  that  differentiation:  her  daughters 
will  remain  in  the  same  compartment.  If  the  cell  had  been  marked 
before  the  creation  of  a  boundary,  her  daughters  would  have  been 
able  to  go  in  any  place  of  the  imaginal  disc,  and  her  descendants 
could  have  been  observed  in  any  of  the  compartments  which  form 
later.  So  a  cell  marked  just  after  the  creation  of  the  anterior 
and  posterior  compartments  will  have  daughters  only  in  the  same 


compartment  as  itself,  but  everywhere  in  this  compartment.  Sub¬ 
sequent  separations  form:  dorsal-ventral,  wing-thorax, . . .  .  At 

each  compartment  formation  something  induces  a  different  commit¬ 
ment  to  the  cells  (and  their  daughters)  according  to  whether  they 
are  in  one  part  or  another  of  the  spatial  domain.  Kauffman  et  al 
[18]  suggested  that  the  cause  could  be  chemical  patterns  of  con¬ 
centrations  in  a  2-species  diffusion-reaction  system: 


Here  ft  is  the  imaginal  disc  and  y  is  a  parameter.  Depending 
on  the  values  of  y,  the  system  (4.1)  possesses  various  stable 
non-uniform  steady  state  (N.U.S.S.)  solutions.  When  such  a 
solution  arises,  those  cells  where  s,  for  example,  is  above 
some  threshold  concentration,  will  differentiate  in  some  way,  the 
others  in  the  alternate  way. 

The  chemical  basis  of  morphogenesis  had  already  been 
suggested  by  Turing  [19]  in  what  is  more  and  more  considered 
as  a  pioneering  paper.  Later  on  Prigogine  and  his  group  [2] 
pointed  out  that  self-organization  could  result  from  the  inter- 
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action  of  diffusion  and  reaction  in  dissipative  structures. 

On  the  other  hand  artificial  enzyme  membranes  (Thomas 
et  al  [8])  are  a  good  tool  to  study,  in  a  well-defined  context, 
basic  phenomena  occuring  in  living  systems:  enzymes  are  catalysts 
of  most  biochemical  reactions,  and  in  living  media  they  are,  in 
general,  immobilized  among  other  protein  molecules  in  structues 
where  diffusion,  in  addition  to  reaction,  play  a  role.  By  link¬ 
ing  enzyme  molecules  to  inactive  protein  molecules,  Thomas  and 
coworkers  prepare  artificial  enzyme  membranes  in  which  it  is 
possible  to  obtain  hysteresis,  oscillations,  pattern  formation 
and  wave  front  propagation.  (Kernevez  et  al,  [20,  21]).  The  fact 
that  pattern  formation  can  occur  in  so  simple  a  context  and  in  a 
structure  which  is  likely  to  occur  in  many  biological  situations 
adds  support  to  Kauffman’s  theory. 

The  so-called  S  -  A  system  models  the  situation  of  an 
active  layer  separated  from  an  external  reservoir  by  an  inactive 
layer  (Figure  4.1). 


The  modelling  assumes  that  in  the  active  layer  the  concentrations 
of  S  and  A  depend  only  on  x  and  y,  coordinates  in  the 
plane  of  ft,  and  that  the  fluxes  of  S  and  A  coming  from 
the  reservoir  to  the  active  layer  through  the  inactive  layer  are 
simply  proportional  to  the  differences  of  concentrations  Sq  -  [S] 
and  Aq  -  [A]  where  Sq  and  Aq  are  concentrations  in  the 
reservoir.  Within  the  active  layer  3  phenomena  occur:  diffusion 
in  the  directions  of  the  plane  of  ft,  with  coefficients  Dg  and 
Da,  S  inhibited  and  A  activated  reaction,  with  rate 

(4.2)  r  =  VM  - L§J -  , 

A  K  +[S](1+^1) 

b  Kss 


and  feeding  from  outside  with  fluxes 


(4.3) 


JS  "  DS 


,  Sr>“  [S] 


and  JA 


,  Aq-[A] 

°aAt 


The  coefficients  in  (4.2)  are,  respectively, 

VM  =  constant,  proportional  to  the  enzyme  concentration 

Ka  =  Michaelis  constant  of  the  enzyme  for  A 

Kg  =  Michaelis  constant  of  the  enzyme  for  S 

KSS  = 


inhibition  constant. 


where  E  is  the  enzyme  uricase,  S  is  uric  acid,  A  is  oxygen, 
and  P  allantoin  plus  other  products. 

The  balance  of  masses  for  S  in  an  elementary  volume  of 
active  layer  gives  (Figure  4.2) 


whence  we  obtain,  by  taking  as  new  units  Kg  for  concentrations. 


L  (i.e.,  some  characteristic  length  of  ft,  for  example,  its 
diameter)  for  length,  9  =  L^/Dg  for  time,  and  dividing  by 
L^dftiKg/e: 


new  coordinates  with  L  as  unit:  x 


Similarly  we  have  for  A 


Finally,  we  have  no- flux  boundary  conditions  on  r  *  3ft 


66 


where  3/3n  is,  say,  the  outward  normal  derivative.  Summing 
up,  the  S  -  A  system  is  governed  by 


(4.8) 


-  As  +  Y [F (s ,a)  -  (sQ-s)]  =  0 
||  -  BAa  +  y (F(s ,a)  -  a(aQ-a)]  =  0 
F(s,a)  =  pas/(l+s+ks2) 


3s 

Tn  " 


0, 


da 


0  • 


We  shall  use  the  notation 


f(s,a)  =  F(s , a)  -  (sQ-s) 
g(s,a)  =  F(s,a)  -  a(aQ-a) 

in  the  following  discussions. 

Let  us  point  out  the  fact  that  all  the  parameters  are 

without  dimension:  Sq  and  Sq,  representing  concentrations  in 

the  environment;  p  which  is  the  ratio  between  2  times;  6 ', 

characterizing  the  access  of  substrates  from  the  outside  to  the 

KA 

active  layer  through  the  inactive  layer;  xr~  »  characteristic 

VM 

time  of  the  enzyme  reaction;  a  and  $»  ratios  of  diffusion 
coefficients;  k,  ratio  of  concentrations;  and  y»  ratio  of 
9  and  0',  characteristic  times  of  diffusion  within  the  active 
layer  and  of  transport  to  the  active  layer. 


In  applications  only  k  and  g  are  imposed:  k  of  the 
order  of  0.1  or  1,  6  of  the  order  of  5.  The  other  parameters, 

s0,a0’  a»  p  an<*  Y»  are  at  our  ^isPosal* 

It  can  be  shown  (Kernevez  [20])  that  these  parameters  can 
be  chosen  in  such  a  way  that  the  so-called  isocline  curves 
f(s,a)  *  0  and  g(s,a)  =  0  (I  igure  4.3)  have  only  one  common 
point  (s,a),  a  stable  steady  state  (node)  of  the  dynamical 
system 

(4.9)  £f  ♦  f(s,a)  =0  +  g(s,a)  =  0 


Equivalently,  the  nature  of  the  steady  state  (5, a)  for  (4.9) 


can  be  expressed  by  saying  that  \1  >  0  and  X2 
being  the  eigenvalues  of  the  linearized  operator 


around  (s,a).  Let  us  calL 


det(O) 


where  F 


(s,a)  is  a  stable  node  iff  tr  (0)  -  4  det(0)  >  0  ,  det  (0)  >  0 
and  tr(0)  >  0. 


We  shall  see  in  the  following  that  the  first  inequality 
will  be  a  consequence  of  others.  Therefore  we  only  require  the 
following  conditions 


det(0)  *  aF  +  F  ♦ 


Now  one  can  observe  that  (s,a)  is  a  trivial  solution  not  only 
of  (4.9)  but  also  of  (4.8).  This  is  due  to  the  fact  that 
s(x,y,t)  s  s,  a(x,y,t)  =  a,  profiles  of  concentration  constant 

3  S 

in  time  and  spatially  uniform,  satisfy  =  0,  *  0,  -As  =  0, 

-Aa  =  0,  f(s,a)  =  0,  g(s,a)  =  0,  =  0  and  *  0. 

As  (s,a)  is  stable  for  (4.9),  one  can  ask  if  (s,a)  is 
still  stable  for  (4.8).  All  the  parameters  being  fixed,  except 
y,  we  shall  see  that  the  answer  is  yes,  if  y  is  small  enough, 
and  no,  if  y  is  large  enough.  (y  small  means  0  small  with 
respect  to  8*,  i.e.,  diffusion  of  S  within  the  active  layer 
fast  in  comparison  with  the  transport  of  S  from  the  outside  to 
the  active  layer).  In  fact,  the  situation  is  more  complex  since, 
as  we  shall  see,  when  y  is  increasing,  (s,a)  can  lose  its 
stability,  then  recover  it,  and  this  can  occur  several  times, 
before  the  point  finally  becomes  and  remains  unstable. 

| 

Finally,  let  us  suppose  that  in  the  following  all  the 

parameters,  except  y,  will  be  kept  fixed,  y  will  be  the  bi- 

furcation  parameter.  Let  us  recall  that  y  ■  ^-0',  so  that  if 

US 

Dg  and  0'  are  kept  fixed,  y  is  proportional  to  meas(ft). 
Therefore,  one  can  view  (4.8)  as  the  modelling  of  an  expanding 
disc  (like  an  imaginal  disc)  after  normalization  of  the  spatial 
domain  to  a  domain  8  with  diameter  1,  the  size  of  the  disc 
being  determined  by  the  coefficient  y.  If  the  size  of  the  spatial 
domain  does  not  vary,  y  may  change  for  other  reasons,  for  ex- 

t 

ample,  because  0'  increases,  and  since  0'  ■  LjL2/Ds,  this  can 

t 

be  due  to  a  decrease  of  Dg. 
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M(U) 

(4.14) 


N(U)  -  F(s+u,a+v)  -  F(s,a)  -  Fcu  -  F  v 

5  cl 

(N(U)  -  0 (u2+v2)) . 


It  can  be  shown  (Meurant  and  Saut  [26])  that  the  operators 
I*Y  ■  Lq  +  yLj  and  -  yM  possess  a  number  of  desirable 
properties  such  as: 
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Lemma  4.1. 
(i) 

(ii) 

(iii) 


(iv) 


For  all  y  6R+ 

Ly  is_  a  closed  linear  operator  with  domain 

0-((u,v)  G  H2(fi)2,  =  0,  -  0}  in  H  -  L2  (n) 2 . 

Ly  is  m-  sectorial  with  half  angle  ir/4. 

The  resolvent  of  is  compact  (the  spectrum  of 

Ly  i^  therefore  discrete  and  only  composed  of 

eigenvalues  having  a  finite  multiplicity) . 

L  is^  the  infinitesimal  generator  of  an  analytic 
-Lyt 

semigroup  e  which  satisfies 


-L  t 

I  le  Y  I  liz^(K,D)  1  ^T77»  t:  e  ]0,T[,  c  >  0, 
where  K  =  H1(fi)2. 

(v)  (Ly)  is.  a  holomorphic  family  of  type  (A) . 


Lemma  4.2. 

(i)  My  K  is^  analytic 

(ii)  I |My(0) | lK  <  C||U||2 

they  are  satisfied 
crosses  the 

remaining  part 
0,  then  we  have: 


The  reason  we  give  these  properties  is  that  when 

and  when  Yn  is  such  that  the  spectrum  of  L_ 
u  Y0 

imaginary  axis  by  the  simple  eigenvalue  0,  the 
of  the  spectrum  lying  in  the  half  space  Re(z)  > 
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Theorem  4.1.  There  exists  in  a  solution  U(y)  *  [u(Y), 
v(Y)]T  %  (0,0]T  of  LyU  +  M  (U)  -  0  which  is  defined  in  a 
neighborhood  of  yq.  Moreover,  if^  the  bifurcation  is  bilateral 
(U(y)  defined  on  the  left  and  right  sides  of  Yg) »  U(Y)  is 
analytic  in  Y,  and  if  the  bifurcation  is  one-sided,  U(y)  is 
analytic  in  Iy-YqI1^.  If^  the  bifurcation  is  bilateral ,  there 
is  "exchange  of  stabilities" .  If  the  bifurcation  is  one-sided 
U(Y)  is  stable  or  not  according  to  whether  it  is  on  the  side 
where  the  trivial  steady  state  is  unstable . 

From  a  practical  viewpoint,  we  must  now  effectively  study 
the  eigenvalues  of  Ly,  to  see  whether  we  are  in  the  case  where 
Theorem  4.1  is  applicable;  this  will  be  the  subject  of  this  para¬ 
graph.  We  then  shall  have  to  elucidate  the  nature  of  the  bifur¬ 
cation  (bilateral  or  one  sided,  on  which  side);  this  will  be  done 
in  §3  by  a  modified  perturbation  analysis. 

Now  the  problem  is:  we  know  that  U  =  0  is  stable  or  not 
according  to  whether  the  spectrum  of  Ly  lies  in  Re(z)  >0  or 
at  least  one  of  its  eigenvalues  lies  in  Re(z)  <0.  So  we  are 
seeking  X  and  $  such  that  (the  eigenvalue  problem) 

(4.15)  Ly4»  -  X*. 


We  can  look  for  $  in  the  form 


where  wn  is  the  n  eigenfunction  of  -L  on  Q  subject  to 
no-flux  boundary  conditions: 


Substituting  (4.16)  into  (4.15),  we  obtain 


This  admits  a  solution  iff  the  determinant  of  the  matrix  in 


(4.18)  is  null,  that  is 


X  -  tr(n)X  +  det(n)  =  0 


Here 


* 

~  1  “ 

(4.16) 

♦n  * 

_  Mn  _ 

Therefore, to  each  eigenvalue  pn  of  -A  (n=0,l,2, . . .),  there 
will  correspond  2  eigenvalues  of  Ly,  X "  and  X*,  solutions 
(4.19) ,  and  to  each  one  will  be  associated  a  value  M~  by 


thus  determining  the  corresponding  eigenvector  4>, 


Proposition  4.1.  Under  the  assumptions 


then  the  solutions  X  of  (4.19)  are  real  numbers  (X“  <  X*) 

11  '  n  n 


Proof 
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But  A^  +  A*  =  tr(n)  =  ( 8+1)  Un+Ytr^0^  >  0  from  C4*11)*  s0  that 
at  least  A*  >  0 .  In  order  for  A~  to  be  negative,  we  must 
have  det(n)  =  A' A*  <  0.  But  det(n)  =  y2T(-^0,  so  that  we 
are  led  to  make  the  following  assumption  on  the  polynomial  T: 

(4.24)  T(z)  admits  2  positive  roots  0  <  z'  <  z". 

y  y  y 

Now  we  have  Xn<0<^z'  <Y?1  <  z"  <£=>  <  Y  <  ~r<F>y  €  In  = 

]z""1vn*  z,-1Vn[.  Thus  we  have  the  following  result. 


Theorem  4.2.  Suppose  that 


(i) 

tr (0)  =  Fs  +  1  +  Fa  + 

a  >  0 

(ii) 

det(0)  =  aFg  +  Fa  +  a 

>  0 

(iii) 

a  >  1,  $  >  1 

(iv) 

(e(Fs+l)  -  (Fa+a))2  + 

4BFaFs  >  0 

(v) 

BCV1)  +  Fa  +  a  <  0. 

Then  there  exists  a  sequence  of  intervals  In  =  ]z"  ,  z' 

such  that 


(1°) 

if 

Y  G  U 
n>l 

then 

(s,a) 

is  unstable 

(2°) 

if 

Y  e  U 

n>l 

T„ 

then 

(s,a) 

is  stable. 

of  L. 


In  a  similar  manner  we  could  define  the  eigenfunctions 
T. 


n 


i 


m 
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(4.25) 


.  T  +  +  + 

Y  n  n  vn 


with  the  same  eigenvalues  XR  but  this  time,  if 


(4.26) 


1±  i"1  ] 

**  =  ;  n*  |Wn 

L  nJ 


then  Nr  is  determined  by  (in  place  of  (4.21)) 


(4.27) 


%  *  Y(FS*1)  -  X-  .  yFaN-  =  0. 


It  can  be  easily  checked  that 


(4.28) 


+  + 
(<f>"t'l'_) 


(4)±.^±) 

vrn  *  Mir 


1  +  M*N*  i  0 
n  n 


n  >  0 


=  0 


=  0 


m  t  n. 


Here,  if  f  =  [fx,f2)T  €  L2(«)2  and  g  =  [g1,g2]T  G  L2(fl)2, 
(f,g)  means 


(4.29) 


(f »g)  -  {  (flgl+f2g2)dn. 
n 


Finally,  we  shall  use  later  the  property  that  the  <J>n  constitute 
a  basis  of  L2(fi)2;  let  fei.2(Jl)2.  Then 
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(4.30) 


f=  l  «>  (f.tf  " 

>u 


Exercise :  Prove  (4.28). 

Hint:  From  (4.21)and  (4.-7)  we  have  FaM“  =  FsN“  so  that 

1  *  MnNn  *  '  1  ♦CM,1*Y(Fs.l)-X^/CB,in*YCFa.,.)-X*) 

•  (tr(n)  -  2X*)/(8VYtFa.«)-X*)  i >  0. 

Similarly  one  can  prove  the  following  which  will  be  useful 

later: 

1  +  3MnNn  =  2eY^  ‘  T?(Fa+0l+6(Fs  +  1))'(e+1)Xn)/(:8VY(Fa+a)'X^)- 
Now  (^n’^n^  =  |wn(1+M^Nn)dJJ  "  1+MnNn  * 

n 

(<P±,'P±)  *  M  w  w  d^)  (1+M±N±)  =  0  since  j  w_w_dfi  =  6 . 

n’Mir  \J  m  n  /  v  n  nr  J  m  n 


'Q 

+  ±  + 


+  + 


+  T  +. 

—  T  4  .1.  Y 


mn 


Finally,  X'($n,*n)  =  (S*n^n]  =  ^n>LA]  =  Xn(*n»*n>  50  that 


<CO  -  °* 


In  §3  the  eigenvalues  X*  and  eigenvectors  4>n  and 

will  be  those  of  Lv  ,  where  Yn  is  a  value  of  the  parameter 

Y0  u 


y  for  which  the  trivial  steady  state  U 


is  losing  stability 


bv  entering  into  In^.  So  that  =  0,  CLq+YqL^) 4>n^  =  0, 


•*n  > 
*0 


1  ♦  M‘  N"  -  tr(nn) / (8u  +Yn(F  +a))  >  0  and  a  quantity 
n0  n0  n0  u  a 


which  shall  be  needed  later 


according  to  whether 


4.3.  Bifurcation  analysis 


We  are  going  to  use  the  modified  perturbation  method  as 


indicated  in  H.B.  Keller  [23] 


Statement  of  the  problem.  Let  yn  be  a  critical  value  of  y 


for  which  the  trivial  steady  state  U  =  0  of  L^U  +  M^fU)  =  0 
is  losing  its  stability  when  y  enters  into  the  interval  I 

n0 

(i.e.,  y_  /yn  =  z'  or  z",  depending  upon  which  extremity  of  I 

0  0  ,  .  .  n0 
Y  crosses).  Let  be  the  eigenvalues,  and  4>n(resp.  ¥“)  be 

T 

the  eigenvectors  of  L  (resp.  L  ).  We  are  seeking  nontrivial 

Y0  Y0 

solutions  of 


LU  +  yM(U)  *  0 


when  y  stays  in  some  neighborhood  of  Yq.  For  that  we  consider 
U  and  y  as  smooth  functions  of  a  parameter  e: 


(4.32) 


U  =*  U(e)  and  y  ■  y(e) 


satisfying  the  "inflated"  problem 


(4.33) 


(L0+y(e)L1)U(e)  +  y(e)M(U(e))  «  0 


(U(e),^  )  *  e(l+M"  N"  ) 
n0  0  0 


and  look  for  Taylor  expansions 


(4.34) 


U(e)  -  eU'(O)  +  £e2U"(0)  +  0(e3) 

Y(e)  -  Y0  +  ey'(O)  +  ^e2y"(0)  +  0(e3) 


First,  we  differentiate  equations  (4.33)  with  respect  to  e 


'  Y'(e)L1U(e)  +  (L0+y(e)L1)U'(e)  ♦  y'(e)M(U(e)) 

(4.35)  J  +  y(e)M’ (U(e))U’(e) 

(u*(e),<»:  )  - 1  +  m;  n:  . 
no  “o  no 


Then,  in  (4.35)  we  put  e  ■  0 


(4.36)  (L0+y0L1)U»(0)  -  0,  (U»  (0)  .f^)  -  1  ♦  M~^N~ 


since 


U(0)  -  0,  M(0)  -  0  and  M' (0)  -  0. 
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From  (4.36)  it  results  that 


(4.37) 


U'(0)  -  4> 


nf 


Now  we  differentiate  equations  (4.35)  with  respect  to  e,  and 
omit  e  everywhere  to  simplify  notation: 


(4.38)  J 


Y'^U  +  2y'L1U’  +  (Lg+yL^U"  ♦  y"M(U)  +  2y'M»(U)U’ 
+  yM"(U)U’2  +  yM'(U)U”  =  0 


(U",*‘  )  -  0  . 
n0 


In  (4.38)  we  again  put  e  *  0: 


(4.39) 


♦  (Lq+YqL^U'^o)  +  y0m,’(°)^o  -  0, 


(U"(0),*n  )  -  0 

n0 


From  the  Fredholm  alternative,  (4.39)  will  have  a  solution 
foT  U"(0)  if  and  only  if 


(2y'  (0)L.4"  +  Y0M"(0)*;2,i|»;  )  -  0. 

1  n0  u  n0  n0 


(4.40) 
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This  relation  gives  y'(0): 


I 

I 


(4.41) 


Y*  CO)  -  -y0 


(M"(0)  ) 

n0  n0 

2(li*;  ) 

i  nQ  nQ 


Recall  that  at  the  end  of  §2  we  did  calculate  the  denominator  of 
(4.41).  On  the  other  hand  the  numerator  must  be  interpreted  in 


the  following  way: 


M(U)  -  N(U) 


[:] 


(F(s+u,a+v)  -  F(s,a)  -  Fgu  -  Fav) 


[:]■ 


If  H  -  [h1,h2]  , 

(M'(U),H)  -  (Fg(s+u,a+v)h1  +  Fa(s+u,a+v)h2  -  F^  - 


M" (U)  (H,K)  -  (f”s  (s+u,a+v)h1k1  +  Fsa(I+u,a+v)  (h1k2  ♦  h^) 


+  Faa(§+u»2+v)h2k2) 


~l“ 

_1_ 


M'  ”  (U)(H,K,L)  -  (Fsss  h1k1i1  +  Fssa  (hjk^j  ♦  h1k#»1  ♦  hjk^) 


+  Fsaa  (hlk242  +  h2kl42  +  h2k24l5  +  Faaa  h2k2425 


From  this,  and  the  fact  that  F  "0,  we  find 


In  the  case  where 


and  the  bifurcation  will  be  unilateral  (Figure  4.4  (a),  (b)). 
Otherwise,  (this  will  be  the  general  case)  the  bifurcation  will 
be  2  sided  (Figure  4.4  (c),(d)). 


So  the  only  case  where  one  needs  further  information  is  the  case 
where  (4.43)  holds.  In  that  case  we  would  like  to  know  on  what 
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side  of  Yg  the  bifurcated  branch  lies,  i.e.,  know  y"(0). 
When  (4.40)holds  U"(0)  is  given  explicitly  by 


(4.44)  U"(0)  - 


■r0  l  -I  ,♦*)  — Vj  ♦* 

#0  \ 


Now  in  order  to  get  U"  (0)  we  differentiate  once  more  in  (4.38) 
and  set  e  »  0  (recall  that  Y'(0)  =  0): 


3y"(0)L14;o  +  (Lg+YgL^U'^CO)  +  YgM"'(0)*;o  +  3YgM"  (0)  (0)  -  0, 


The  condition  for  existence  of  U,,'(0)  is 


Yg  .3 

w  •  t  •  /  AY  a 


(4.45)  Y  (0)  -  -  - .  ■.  (CM’ “(0)4  ,*  )  +  3(M"(0)4  U"(0)^_  )) 

3(Ll^n  »♦«  )  0  0  0  0 

i  nQ  nQ 


and  this  quantity  can  be  calculated. 


Remark .  In  the  case  of  a  2-  or  3-  dimensional  domain  ft  with 

arbitrary  shape,  there  is  no  reason  for  (4.43)  to  hold.  However, 

f1  3 

in  the  1-dimensional  case,  where  w_(x)  -  cos  nirx,  w  dx  -  0. 

“  Jn 


'dx  •  0. 


Also  when  ft  has  1  axis  of  symmetry,  say  0  -  y  suppose  that  w 
is  an  eigenfunction: 


■Aw(x,y)  -  Xw(x,y)  ,  j  w' 


dft  -  1 


.  •  .  .  "j  m  -rf-dKafa  -||,  -iy-  - 1 


and  consider  v(x,y)  *  w(-x,y). 


-Av(x,y)  -  -Aw(-x,y)  -  Xw(-x,y)  -  Xv(x,y),  J  v2dft  -  1 

ft 

Then  v(x,y)  -  ±w(x,y),  i.e.,  w(-x,y)  -  ±w(x,y),  so  that  there 

is  the  possiblity  for  w  to  satisfy  w(-x,y)  *  -w(x,y).  In 

that  case  j  w3dft  *  0  and  (4.43)  holds, 
ft 

The  case  ft  »  ]0,1[  is  not  unrealistic.  It  corresponds 
to  the  modelling  of  an  ellipsoid  with  large  axis  ratio  when  one 
is  interested  in  the  first  eigenmodes.  One  can  explain  the 
differentiation  of  regions  in  an  egg  by  a  steady  pattern  as 
pictured  in  Figure  4.5  (Kauffman  et  al  [18]).  See  also  §5. 


Figure  4.5 


4.4.  Numerical  analysis 


Two  kinds  of  methods  are  at  our  disposal: 

(i)  methods  in  which  the  steady  states  are  obtained 
by  solving  the  evolution  equations  until  one 
estimates  that  the  system  has  attained  a  (stable) 
steady  state;  and 

(ii)  methods  following  branches  of  (stable  or  unstable) 
steady  state  solutions. 

The  second  type  of  method  was  discussed  in  Chapter  2. 


Solution  of  the  evolution  equations 


Discretization  with  respect  to  time 


Discretization  with  respect  to  space 


By  finite  differences  in  the  1-dimensional  case 
(by  finite  elements  in  the  2-dimensional  case 
(Figure  4.6))  (4.46)  can  be  written 


•+ya)M)  A1 
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Figure  4.6 


where 

K  ■  stiffness  matrix  (corresponding  to  -6) 

M  -  mass  matrix  (corresponding  to  I) 
fn  -  vector  -  y(F(sn,an)  -  sQ) 
gn  -  £  -  y(F (sn,an)  -  oa0). 

For  the  practical  formulation  and  implementation  of  the 
finite  element  method  we  refer  to  Bathe  and  Wilson  [15] .  Let  us 
give  more  details  on  the  calculation  of  the  "stiffness"  and  "mass" 
matrices  K  and  M  in  the  case  of  a  surface  (a  planar  domain 
being  a  particular  case).  The  surface  is  divided  into  curved 
8-nodes  quadrilaterals  (see  Figure  4.7).  The  nodes  are  defined 
by  their  coordinates  for  node  number  i  of  the 

quadrilateral,  1  <  i  <  8).  We  call  H^(r,s)  the  "shape  functions" 


87 

H1(r,s)  =  j(l+r)  (1+s)  -  ^H5(r,s)  -  £«8(r,s) 

H2(r,s)  85  *(l-r)(l+s)  -  X(r's)  "  7^5 Cr,s) 

Hj(r,s)  -  J(l-r)(l-s)  -  ^y(r,s)  -  ^«6(r,s) 

H4(r,s)  -  j(l+r)  (1-s)  -  ^Hg(r,s)  '  iM**5) 

H5(r,s)  -  y(l-r2)  (1+s) 

H6(r,s)  -  £(l-s2)(l-r) 

H7(r,s)  -  £(l-r2)(l-s) 

H8(r,s)  *  yCl-s2)  (1+r) . 


Figure  4.7 
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The  part  of  the  surface  delineated  by  the  quadrilateral  is  approxi¬ 
mated  by  the  surface  whose  parametric  representation  is 


x 


(4.47a)  < 


8 

l  x-H. (r,s),  y 
i=l  1  1 


-1  <  r,  s  <  +1. 


8 

I  y^H. (r,s)t  z 
i-1  1  1 


8 

l  ziHi(r,s) 
i  =  l  1  1 


The  elementary  mass  and  stiffness  matrices  are 


(a) 


B(m)  D(m)B(m)ds(m)  and  M(m)  f  H(^dSW 

Js(m)  J  (m) 


s(m) 

where  dS<m),B(m), 


i  Cm) 


and  are  defined  below. 

To  each  pair  (r,s)  6  [-1,+1]  *  [-!,♦!]  corresponds  a 


point  (x,y,z)  on  S 


(m) 


by  equations  (4.47a)  and  we  call 


Hi"^(x,y,z)  =  H^r.s)  1  <  i  <  8. 

Thus  is  a  function  of  the  point  (x,y,z)  on  the  surface 

S(m).  denotes  the  vector 


H 


(a)  _ 


[H^m) ,H^m) 


.H«] 


Similarly  we  define 


and  are  8xl  and  8x2  matrices  defined  on 

Now,  following  equations  (1.33)-(1.36)  of  Chapter  I,  we  define 


3H .  3H. 


dx2  *  ay2  *  a,2  -  (.1  *♦£*)) 

1=1  ' 


3H .  3H . 


3H.  3H. 


2  2 
alldr  +  2a^2drds  +  a22ds 


all  =  V1  »  a12  V1V2 *  a22  "  V2 

♦  f  8  3H.  3H .  3H,T 

vi  s  [J^iTT*  ly±  TT>  lz i  Trj 

T  3H.  3H,  3H.1T 

V2  "  fxi  7s~*  Eyi  7s-’  ^z.i  7s-] 


Then  we  define 


a  "  alla22  "  a12 


a11  =  a22/a  ,  a22  =  an/a  *  ®12  *  "ai2^a 
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dS(ra) 

D(m) 

K(m) 

mW 


a1/2drds 


> 


1 

22  I  » 


;-l<r,s<+i 
^  ~ 1<T , S <+  1 


B(m)  jjCmig  tra)a^^2drds , 

H(m)T  HW.i^drds. 


Using  this  method  and  the  program  SSPACE  (subspace  iteration 
method)  of  Bathe  and  Wilson  [15],  the  first  eigenvalues  and  eigen¬ 
functions  of  the  Laplace  Beltrami  operator  -Ag  on  the  surface 
of  an  elongated  ellipsoid  with  the  form  of  an  egg  can  be  obtained 
(Bunow  et.al.  [24]);  (See  Figures  4.8,  4.9). 
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(iii)  separation  between  left  and  right  sides 

etc « .  • 


Surfaces  where  the  diffusion  coefficients  are  varying 
(because  of  a  layer  of  cells  of  variable  size  for  instance)  model 
■ore  closely  the  shape  of  the  real  imaginal  discs  (Figure  4.10). 
It  is  thus  possible  to  solve  the  S-A  system  equations  on  such 
surfaces  and  study  the  sequence  of  stable  patterns. 


Section  of  an  imaginal  disc 


states  (N.U.S.S.)  froa  the  nuaerical  results  in  the  1 -dimensional 
case,  where  one  has  situations  as  in  Figure  4.12. 
h  A  I  V  3- structuration 


4 -structuration 


-  stable 

-  unstable 

Figure  4.12 

Hysteresis  cycle  when  y  first  increases,  then  decreases,  the 
steady  state  following  the  path  1-2-3-4-S-6-7 . 


Reaction  -  diffusion  model  of  Murray  et  al  [25],  for  the  hydranth 
regeneration  in  the  marine  hydroid  Tubularia. 

Pigment  distribution  presages  hydranth  regeneration  of  an 
amputated  hydranth  in  the  marine  hydroid  Tubularia.  Murray,  et  al, 
[25],  suggest  that  such  a  distribution  could  result  from  a  reaction 
diffusion  system.  In  the  marine  hydroid  Tubularia  regeneration  of 
an  amputated  hydranth  takes  place  in  the  order  of  one  and  a  half 
days.  The  sequence  of  clearly  observable  states  is  illustrated 


schematically  in  Figure  4.13.  In  Figure  4.13  (a),  (b)  and  (c) 
the  shading  and  serrated  lines  represent  clearly  seen  reddish 
coloured  pigment.  The  regular  periodic  configuration  in  Figure 
4.13(c)  presages  the  tentacle  distribution  in  Figure  4.13  (d) . 
Murray  considers  the  body  of  the  Tubularia  to  be  represented  by 
a  circular  cylinder  of  internal  radius  r^  and  external  radius 
r,  and  considers  the  reaction- diffusion  mechanism  model  to  be 


similar  to  the  s  -  a  system  with  a  =  1.5,  8  *  100,  k  *  0.1, 

Sq  -  92,  Sq  ■  65,  p  »  18.5  (  s  =  10,  a  ■  9.4)  and  various  y. 
Either  s  or  a  is  related  to  the  pigmentation  and  their 
spatial  structure  indicates  the  colouration. 

The  axial  distribution  in  Figure  4.13  (a)  and  (b)  is 
explained  by  the  1-dimensional  version  of  the  s  -  a  system. 

The  circumferential  distribution  in  Figure  4.13  (c)  can 
be  produced  by  a  2-dimensional  version  of  the  s  -  a  system  in 
the  geometry  obtained  by  a  cross  section  of  the  finite  wall 
cylinder. 


patterns  in  a  cross  section  of  the  finite  wall  cylinder 
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The  typical  length  L  is  .  The  reaction-diffusion  domain  is 
then  given  by  1  <  r  <  6  (fi-^/rj)  and  the  equations  are 


(4.48) 


3s 

Jt 


3a 

7F 


-  ♦  7  jf  *  7  *  Yffs-a)  '  0 

•  6C0  *  ?  H  *  7  &  *  T‘Cs’,)  ■  0 


with  the  boundary  conditions 


(4.49) 


°  for  r  -  1.  r  -  t. 


Murray  notes  that  the  number  of  tentacles  is  independent  of  over¬ 
all  size,  if,  all  other  parameters  being  constant,  the  ratio 
r2 

~  $  is  kept  fixed.  This  could  explain  the  regulating  feature 

of  Tubularia  that  approximately  the  same  number  of  tentacles  appear, 
roughly  12-18,  irrespective  of  size. 
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